#include "stdafx.h"
#include "Task.h"
#define _USE_MATH_DEFINES
#include <cmath>

Task::Task( TaskType _t, int _x, int _y, double acc, int steps ): type(_t), size_x(_x), size_y(_y), accuracy(acc), num_steps(steps)
{
	double x,y;
	matrix = new double[(size_x + 1)*(size_y + 1)];
	diff = new double[(size_x+1)*(size_y+1)];
	solve(size_x, size_y, matrix);	

	if (type == TaskTypeMain)
	{
		matrix2 = new double[4*(size_x+1)*(size_y+1)];
		solve(size_x*2, size_y*2, matrix2);
		for(int j=0; j<size_y+1; j++)
			for (int i=0; i<size_x+1; i++)
				diff[(size_x + 1)*j + i] = fabs(matrix[(size_x + 1)*j + i] - matrix2[2*(2*size_x + 1)*j + 2*i]);
	}
	else
	{	
		matrix2 = new double[(size_x+1)*(size_y+1)];
		for(int j=0; j<size_y+1; j++)
			for (int i=0; i<size_x+1; i++)
			{
				x = (3.0/size_x)*i; 
				y = ((double) 1.0*j)/((double) size_y);
				matrix2[(size_x + 1)*j + i] = sin((double) x*y*y)*sin((double) x*y*y);
			}
		for (int i=0; i<(size_x+1)*(size_y+1); i++)
			diff[i] = fabs(matrix[i] - matrix2[i]);
	}
}


void Task::solve(int size_x, int size_y, double* matrix)
{
	double* f = new double[(size_x + 1)*(size_y + 1)];
	double x,y;
	bool flag = false;
	double h2 = size_x*size_x/9.0;
	double k2 = size_y*size_y;
	double a2 = 2 * (h2 + k2);

	double eps_max = 0;
	double eps_cur = 0;
	int steps = 0;


	// Initilize 

	if (type == TaskTypeMain)
	{
		for (int i=0; i<size_x+1; i++)
		{
			x = (3.0/size_x)*i; 
			matrix[i] = cosh(x*x - 3*x) - 1;
			matrix[size_y*(size_x + 1) + i] = 0;
		}
		for (int i=0; i<size_y; i++)
		{
			y = (1.0/size_y)*i;
			matrix[(size_x + 1)*i] = sin(M_PI*y)*sin(M_PI*y);
			matrix[(size_x + 1)*i + size_x] = 0;
		}
		for(int j=0; j<size_y+1; j++)
			for (int i=0; i<size_x+1; i++)
			{
				x = (3.0/size_x)*i; 
				y = (1.0/size_y)*j;
				f[(size_x + 1)*j + i] = (-1)*cosh(x - y);
			}
	}
	else
	{
		for (int i=0; i<size_x+1; i++)
		{
			x = 3.0*i/size_x; 
			matrix[i] = 0;
			matrix[size_y*(size_x + 1) + i] = sin(x)*sin(x);
		}
		for (int i=0; i<size_y+1; i++)
		{
			y = 1.0*i/size_y;
			matrix[(size_x + 1)*i] = 0;
			matrix[(size_x + 1)*i + size_y] = sin(3.0*y*y)*sin(3.0*y*y);
		}
		for(int j=0; j<size_y+1; j++)
			for (int i=0; i<size_x+1; i++)
			{
				x = (3.0/size_x)*i; 
				y = (1.0/size_y)*j;
				f[(size_x + 1)*j + i] = 2.0*x*sin(2.0*x*y*y) + 2.0*y*y*y*y*cos(2.0*x*y*y) + 8.0*x*x*y*y*cos(2.0*x*y*y);
			}
	}

	for (int j=1; j<size_y; j++ )
		for (int i=1; i<size_x; i++)
			matrix[(size_x + 1)*j + i] = 0;

	while (!flag)
	{
		eps_max = 0;
		for (int j=1; j<size_y; j++ )
			for (int i=1; i<size_x; i++)
			{
				double v_old = matrix[(size_x + 1)*j + i];
				double v_new = k2*(matrix[(size_x + 1)*(j+1) + i] + matrix[(size_x + 1)*(j-1) + i]) + h2*(matrix[(size_x + 1)*j + i+1] + matrix[(size_x + 1)*j + i-1]) - f[(size_x + 1)*j + i];
				v_new = v_new / a2;
				eps_cur = fabs(v_old - v_new);
				if (eps_cur > eps_max)
					eps_max = eps_cur;
				matrix[(size_x + 1)*j + i] = v_new;
			}
			steps++;
			if ((eps_max < accuracy) || (steps >= num_steps)) 
				flag = true;
	}
	accuracy = eps_max;
	num_steps = steps;
}

double* Task::getSolution()
{
	return matrix;
}

double* Task::getSolution2()
{
	return matrix2;
}

double* Task::getDifference()
{
	return diff;
}

void Task::getReport( double& x_max, double& y_max, double& max, double& acc, int& st )
{
	acc = accuracy;
	st = num_steps;
	double x,y;
	max = 0;
	for(int j=0; j<size_y+1; j++)
		for (int i=0; i<size_x+1; i++)
		{
			x = (3.0/size_x)*i; 
			y = (1.0/size_y)*j;
			if (diff[(size_x + 1)*j + i] > max)
			{
				max = diff[(size_x + 1)*j + i];
				x_max = x;
				y_max = y;
			}
		}
}

Task::~Task()
{
	delete[] matrix;
	delete[] matrix2;
	delete[] diff;
}
